Background

Data are publicly available through Dryad for “Stronger social bonds do not always predict greater longevity in a gregarious primate” publication (Thompson and Cords, 2018). For this publication, Thompson and Cords conducted survival analyses to understand the association between measures of sociability and mortality among blue monkeys in Kenya. Using these data, the analysis described here will attempt to answer two follow-up research questions which expand on the original results presented by Thompson and Cords.

Introduction

In this investigation, blue monkey groups were observed for 9 years to understand the impact of sociability on survival. Among humans and non-human species that demonstrate strong social ties, bonding can be an important factor for longevity, offspring survival, and reproductive rate. Together, these are indicators of fitness.

Social measures collected by Thompson and Cords that are available in their publicly available dataset include bond strength and partner consistency. These were evaluated among each female dyad on an annual basis. For each subject, the bond strength and partner consistency scores for her top three and top six closest partners were averaged. Other available demographic and environmental factors include each subject’s start age and age at censorship, mortality outcome, age at first reproduction, dominance rank, and the number of female groupmates.

For time-varying data, subjects’ social and environmental factors were annually for each study year during which they were observed:

Subject Age at year start Age at year end or censorship Death Bond strength - 3 closest partners Bond strength - 6 closest partners Dominance rank Number of groupmates Consistency - 3 closest partners Consistency - 6 closest partners
adel 4644 5009 0 3.917304 2.513434 0.8666667 16 0.3333333 0.5000000
adel 5010 5374 0 3.603165 2.599880 0.6250000 9 0.6666667 0.8333333
adel 5383 5739 0 6.623989 3.686815 0.6666667 10 0.3333333 0.5000000
adel 5740 5801 1 2.207133 1.241564 0.8181818 10 0.6666667 0.6666667
agne 3673 4038 0 9.184368 5.301906 0.6000000 16 0.6666667 0.5000000
agne 4039 4403 0 2.217799 1.353885 0.2500000 9 0.3333333 0.6666667

In their analysis, Thompson and Cords categorized individuals into one of four groups based on both their bond strength and partner consistency: +/+ (bond strength and consistency above the mean), +/- (bond strength above the mean, consistency below the mean), -/+, and -/-. In both same-year and multi-year models, belonging to the +/- group was found to be the most hazardous for mortality (multi-year HR [referent = +/+], 95% CI: 20.3, 3.4-123). The authors hypothesize that these findings indicate a costliness for female blue monkeys to be part of a strongly-bonded partnership. Females that become highly invested in an inconsistent dyad are at a disadvantage.

Research question

Does a decrease in partnership consistency over the previous year affect the likelihood of mortality for female blue monkeys?

In light of Thompson and Cords’ findings, do individuals who experience a change in closest partnerships, as indicated by a decrease in consistency score from lagges years, experience increased mortality risk as compared to those whose partnership consistency stays the same or increases.

Methods

Simple approach

Here the time-varying dataset will be used. A new variable will be created which indicates a decrease in partnership consistency from the previous year versus partnership consistency that stayed the same or increased. This new indicator variable will be added to the model fit by Thompson and Cords, which included each subject’s social category (+/+, +/-, -/+, or -/-), the number of female groupmates, and their dominance rank during each year of observation.

coxph_mod1 <- coxph(Surv(start.age, stop.age, death) ~ as.factor(soc.cat) + no.groupmates + rank + 
                      as.factor(cons.decrease),
                    data=annual_df)

Time-varying IPW approach

The association between potential confounders and the exposure of interest (loss in partner consistency) will be evaluated. Those that are imbalanced across exposure groups, as well as age at study entry as a time-invariant covariate will be used to predict propensity scores of rloss in partner consistency. Cumulative product inverse probability weights will be calculated for each subject. If necessary, a numerator model will be fit to stabilize the IPWs. A Cox proportional hazards model will be fit with any covariates included in the numerator model, as well as the consistency loss exposure variable. This model will include the cumulative product IPW as the weights option and subject ID for a cluster effect.

Results

Data exploration

## There were  20  deaths among the cohort over the course of this study.

Simple approach

A new variable was be created to indicate a decrease in partnership consistency from the previous year versus partnership consistency that stayed the same or increased. This new indicator variable was be added to the model fit by Thompson and Cords.

coxph_mod1 <- coxph(Surv(time, time2, death) ~ 
                      a.st.co3 + af.groupmates + rank + 
                      entry.years + cons.p3_loss,
                    data=annual_df_ps)
## # A tibble: 7 x 4
##   term             hr low_hr high_hr
##   <chr>         <dbl>  <dbl>   <dbl>
## 1 a.st.co31     1.25   0.209    7.50
## 2 a.st.co32     2.18   0.515    9.27
## 3 a.st.co33     1.96   0.320   12.0 
## 4 af.groupmates 0.965  0.864    1.08
## 5 rank          1.27   0.263    6.12
## 6 entry.years   1.19   1.06     1.34
## 7 cons.p3_loss1 1.54   0.421    5.64
##                 chisq df    p
## a.st.co3      2.62183  3 0.45
## af.groupmates 0.00601  1 0.94
## rank          0.01533  1 0.90
## entry.years   1.22770  1 0.27
## cons.p3_loss  0.11903  1 0.73
## GLOBAL        4.03297  7 0.78

Time-varying IPW approach

Let’s start by checking whether number of groupmates and dominance ranke are, in fact, associated with exposures (partnership consistency and bond strength).

Are the confounders associated with the exposures?

## # A tibble: 4 x 2
##   name             value
##   <chr>            <dbl>
## 1 mean_rank_loss   0.573
## 2 mean_rank_noloss 0.487
## 3 sd_rank          0.323
## 4 stand_diff_rank  0.268

Since the standardized mean difference is >0.1, there is an imbalance in rank across the exposure groups (loss vs. no loss in partner consistency).

## # A tibble: 4 x 2
##   name                     value
##   <chr>                    <dbl>
## 1 mean_groupmates_loss   14.2   
## 2 mean_groupmates_noloss 13.8   
## 3 sd_groupmates           4.79  
## 4 stand_diff_groupmates   0.0799

So number of groupmates is not imbalanced across exposure groups.

## # A tibble: 2 x 9
##   cons.p3_loss n_negneg n_posneg n_negpos n_pospos perc_negneg perc_posneg
##   <fct>           <int>    <int>    <int>    <int>       <dbl>       <dbl>
## 1 0                  36       70       16       64       0.194       0.376
## 2 1                  49        9       24        3       0.576       0.106
## # ... with 2 more variables: perc_negpos <dbl>, perc_pospos <dbl>

The social categories are very imbalanced across exposure groups.

Estimating propensity scores

model_ps <- glm(cons.p3_loss ~ ns(time,df=2) + entry.years + rank_l1 + af.groupmates_l1 +
                  a.st.co3_l1, 
                family = "binomial", data = annual_df_ps)
## # A tibble: 9 x 5
##   term              estimate std.error statistic      p.value
##   <chr>                <dbl>     <dbl>     <dbl>        <dbl>
## 1 (Intercept)        -1.42      0.925     -1.53  0.126       
## 2 ns(time, df = 2)1   2.18      1.08       2.03  0.0428      
## 3 ns(time, df = 2)2  -0.502     0.556     -0.904 0.366       
## 4 entry.years        -0.0648    0.0459    -1.41  0.158       
## 5 rank_l1             0.331     0.514      0.643 0.520       
## 6 af.groupmates_l1    0.0764    0.0387     1.98  0.0482      
## 7 a.st.co3_l11       -2.77      0.511     -5.42  0.0000000611
## 8 a.st.co3_l12        0.759     0.383      1.98  0.0474      
## 9 a.st.co3_l13       -2.73      0.653     -4.18  0.0000288

Calculating IPWs

## # A tibble: 1 x 4
##   `Mean cumulative p~ `Minimum cumulativ~ `Maximum cumulativ~ `Sum cumulative p~
##                 <dbl>               <dbl>               <dbl>              <dbl>
## 1                6.44                1.02                74.0              1744.

Fitting numerator model and calculate stabilized IPW

mod_IPWnum_2 <- glm(cons.p3_loss ~ ns(time, df = 2) + entry.years +a.st.co3_l1, 
                    family = "binomial", data = annual_df_ps)
## # A tibble: 7 x 5
##   term              estimate std.error statistic     p.value
##   <chr>                <dbl>     <dbl>     <dbl>       <dbl>
## 1 (Intercept)        -0.114     0.653     -0.175 0.861      
## 2 ns(time, df = 2)1   2.00      1.06       1.88  0.0595     
## 3 ns(time, df = 2)2  -0.382     0.547     -0.698 0.485      
## 4 entry.years        -0.0727    0.0452    -1.61  0.107      
## 5 a.st.co3_l11       -2.57      0.493     -5.22  0.000000183
## 6 a.st.co3_l12        0.825     0.377      2.19  0.0285     
## 7 a.st.co3_l13       -2.71      0.648     -4.17  0.0000299
## # A tibble: 1 x 4
##   `Mean stabilized IPW` `Minimum stabilize~ `Maximum stabiliz~ `Sum stabilized ~
##                   <dbl>               <dbl>              <dbl>             <dbl>
## 1                 0.992               0.502               1.71              269.

Cox proportional hazards model

coxph_modIPW_tv <- coxph(Surv(time,time2,death) ~
                           cons.p3_loss + entry.years + a.st.co3_l1,
                         weights=w_new_2, cluster=subj,
                         data=annual_df_ps)
## # A tibble: 5 x 6
##   term          estimate std.error robust.se statistic   p.value
##   <chr>            <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
## 1 cons.p3_loss1    0.629    0.651     0.642      0.981 0.327    
## 2 entry.years      0.191    0.0545    0.0477     4.00  0.0000628
## 3 a.st.co3_l11     0.787    0.757     0.654      1.20  0.229    
## 4 a.st.co3_l12     0.600    0.726     0.691      0.868 0.385    
## 5 a.st.co3_l13     0.779    0.843     0.721      1.08  0.280
## # A tibble: 5 x 4
##   term             hr low_hr high_hr
##   <chr>         <dbl>  <dbl>   <dbl>
## 1 cons.p3_loss1  1.88  0.534    6.60
## 2 entry.years    1.21  1.10     1.33
## 3 a.st.co3_l11   2.20  0.609    7.92
## 4 a.st.co3_l12   1.82  0.470    7.06
## 5 a.st.co3_l13   2.18  0.531    8.95
##              chisq df    p
## cons.p3_loss 0.488  1 0.48
## entry.years  1.254  1 0.26
## a.st.co3_l1  1.111  3 0.77
## GLOBAL       4.006  5 0.55

“The correct interpretation of this HR is as the effect of exposure in the exposed. In other words it is the HR comparing what would have happened if the exposed were exposed (what actually happened) to what would have happened if the exposed were unexposed. By contrast remember that the interpretation of the HR in the case of IPW was the HR comparing what would have happened if everyone was exposed to what would have happened if nobody was exposed.”